package scu.maqiang.homogenization;

import scu.maqiang.fes.FES3H81;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H201;
import scu.maqiang.mesh.Mesh3H20;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PossionSOTSProblem3DH {

		public static void main(String[] args) {
			long begin = System.nanoTime();
		    int N = 32;
		    double eps= 1.0 / N;
		    String dir = "PossionSOTS3DH\\";
		    Mesh3H8 cellMesh = new Mesh3H8().cube(32, 32, 32);
		    ScalarFunc regionCellFunc = (xyz, label, param) -> {
	            return (Math.abs(xyz[0] - 0.5) < 0.25 && Math.abs(xyz[1] - 0.5) < 0.25) && Math.abs(xyz[2] - 0.5) < 0.25? 1.0:-1.0;
	        };
		    cellMesh.setDomainLabel(regionCellFunc, null, 10);
		    FES3H81 fsCell = new FES3H81(cellMesh);
		    double[][] FOCS = new double[3][fsCell.GetNdof()];
			double[][] SOCS = new double[9][fsCell.GetNdof()];
			double[][] D = new double[3][3];
			SOTS3DH8201 ss = new SOTS3DH8201(fsCell, null);
			ScalarFunc cf = (xyz, label, param) -> label == 10? 0.01: 1;
//			
			ScalarFunc[] gcf = {cf};
			ss.KxxCellSolution(gcf, null, FOCS, D, SOCS);
//		    
			System.out.println("min: " + MVO.min(FOCS[0]) + "\tmax: " + MVO.max(FOCS[0]));
			System.out.println("min: " + MVO.min(FOCS[1]) + "\tmax: " + MVO.max(FOCS[1]));
			System.out.println("min: " + MVO.min(FOCS[2]) + "\tmax: " + MVO.max(FOCS[2]));
			System.out.println("Homogenized Kxx: ");
			System.out.println(MVO.toString(D));
			System.out.println("min: " + MVO.min(SOCS[0]) + "\tmax: " + MVO.max(SOCS[0]));
			System.out.println("min: " + MVO.min(SOCS[1]) + "\tmax: " + MVO.max(SOCS[1]));
			System.out.println("min: " + MVO.min(SOCS[2]) + "\tmax: " + MVO.max(SOCS[2]));
			System.out.println("min: " + MVO.min(SOCS[3]) + "\tmax: " + MVO.max(SOCS[3]));
			System.out.println("min: " + MVO.min(SOCS[4]) + "\tmax: " + MVO.max(SOCS[4]));
			System.out.println("min: " + MVO.min(SOCS[5]) + "\tmax: " + MVO.max(SOCS[5]));
			System.out.println("min: " + MVO.min(SOCS[6]) + "\tmax: " + MVO.max(SOCS[6]));
			System.out.println("min: " + MVO.min(SOCS[7]) + "\tmax: " + MVO.max(SOCS[7]));
			System.out.println("min: " + MVO.min(SOCS[8]) + "\tmax: " + MVO.max(SOCS[8]));
			cellMesh.toTecplot(dir + "CellFOCSH.dat", FOCS);
			cellMesh.toTecplot(dir + "CellSOCSH.dat", SOCS);

			Mesh3H20 homoMesh = new Mesh3H20(cellMesh);
			FES3H201 fsHomo = new FES3H201(homoMesh);
			ScalarFunc homoCoef = (xy, label, param) -> D[0][0];
			SRMatrix homoA = new SRMatrix(fsHomo.GetNdof());
			ScalarFunc[] MyHomoCoef = {homoCoef};
			fsHomo.assembleStiff(MyHomoCoef, null, BVPType.COMMON, homoA);
			double[] homoRHS = new double[fsHomo.GetNdof()];
			fsHomo.assembleSource(new double[] {10.0}, BVPType.COMMON, homoRHS);
			fsHomo.applyBC_MBN(homoA, Direct.X, 1, 2, 3, 4, 5, 6);
			fsHomo.applyBC_MBN(homoRHS, Direct.X, 0.0, 1, 2, 3, 4, 5, 6);
			NewIterSSolver homoSolver = new NewIterSSolver(homoA, true);
			double[] homoSol = new double[fsHomo.GetNdof()];
			homoSolver.PCGSSOR(homoRHS, homoSol, 1.5, 1);
////			
			//homoA.clear();
			//homoA = null;
			
			Mesh3H8 fineMesh = new Mesh3H8().cube(256, 256, 256);
			ScalarFunc regionFineFunc = (xy, label, param) -> {
	            double xCell = xy[0] * N - (int)(xy[0] * N);
	            double yCell = xy[1] * N - (int)(xy[1] * N);
	            double zCell = xy[2] * N - (int)(xy[2] * N);
	            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25 && Math.abs(zCell - 0.5) < 0.25)? 1.0:-1.0;
	        };
	        fineMesh.setDomainLabel(regionFineFunc, null, 10);
	        FES3H81 fsFine = new FES3H81(fineMesh);
	        SRMatrix fineA = new SRMatrix(fsFine.GetNdof());
	        fsFine.assembleStiff(new ScalarFunc[] {cf}, null, BVPType.COMMON, fineA);
	        double[] fineRHS = new double[fsFine.GetNdof()];
	        fsFine.assembleSource(new double[] {10.0}, BVPType.COMMON, fineRHS);
	        fsFine.applyBC_MBN(fineA, Direct.X, 1, 2, 3, 4, 5, 6);
	        fsFine.applyBC_MBN(fineRHS, Direct.X, 0.0, 1, 2, 3, 4, 5, 6);
	        double[] fineSol = new double[fsFine.GetNdof()];
	        NewIterSSolver fineSolver = new NewIterSSolver(fineA);
	        fineSolver.PCGSSOR(fineRHS, fineSol, 1.5, 1);
//	        
////	    fineMesh.toTecplot("FineSolutionH8.dat", fineSol);
	        
//	        double[] x = MVO.linspace(0, 1, 1001);
//	        double[][] xyz = new double[1001][3];
//	        double[][] value = new double[1001][4];
//	        for(int i = 0; i < 1001; i++) {
//	        	Arrays.fill(xyz[i], x[i]);
//	        }
//	        ss.assembleHeatSOTSSolution(fsFine, fineSol, fsHomo, homoSol, fsCell, FOCS, SOCS, N, eps, xyz, value);
//	        Tecplot.LineXY("SOTSSolutionLine.dat", x, value);

	        double[] HomoFineSol = new double[fsFine.GetNdof()];
	        double[] FirstSol = new double[fsFine.GetNdof()];
	        double[] SecondSol = new double[fsFine.GetNdof()];
//	          
	        ss.assembleHeatSOTSSolution_Cube(fsFine, HomoFineSol, FirstSol, SecondSol, fsHomo, homoSol, fsCell, FOCS, SOCS, N, eps);
//	        
	        System.out.println("Homo Solution");
	        System.out.println("Homo min: " + MVO.min(HomoFineSol) + "\tmax: " + MVO.max(HomoFineSol));
	        System.out.println("First Order Approximation: ");
	        System.out.println("FOA min: " + MVO.min(FirstSol) + "\tmax: " + MVO.max(FirstSol));
	        System.out.println("Second Order Approximation: ");
	        System.out.println("SOA min: " + MVO.min(SecondSol) + "\tmax: " + MVO.max(SecondSol));
//	        
	        double[][] wholeSol = new double[4][];
	        wholeSol[0] = HomoFineSol;
	        wholeSol[1] = FirstSol;
	        wholeSol[2] = SecondSol;
	        wholeSol[3] = fineSol;
	        fineMesh.toTecplot(dir + "SOTSSolutions.dat", wholeSol);
	        double L2Norm = fsFine.computeL2Norm(fineSol);
	        double[] error0 = MVO.add(fineSol, -1.0, HomoFineSol);
			double[] error1 = MVO.add(fineSol, -1.0, FirstSol);
			double[] error2 = MVO.add(fineSol, -1.0, SecondSol);
			System.out.println("L2RError:");
			System.out.println("Error 0: " + fsFine.computeL2Norm(error0) / L2Norm);
			System.out.println("Error 1: " + fsFine.computeL2Norm(error1) / L2Norm);
			System.out.println("Error 2: " + fsFine.computeL2Norm(error2) / L2Norm);
//			
			double H1SemiNorm = fsFine.computeH1SemiNorm(fineSol);
			System.out.println("H1SemiRError:");
			System.out.println("Error 0: " + fsFine.computeH1SemiNorm(error0) / H1SemiNorm);
			System.out.println("Error 1: " + fsFine.computeH1SemiNorm(error1) / H1SemiNorm);
			System.out.println("Error 2: " + fsFine.computeH1SemiNorm(error2) / H1SemiNorm);		
			
			System.out.println("LInftyRError");
			System.out.println("Error 0: " + MVO.maxA(error0) / MVO.maxA(fineSol));
			System.out.println("Error 1: " + MVO.maxA(error1) / MVO.maxA(fineSol));
			System.out.println("Error 2: " + MVO.maxA(error2) / MVO.maxA(fineSol));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
		}

}
